# Code for Table J1: Main Effects

rm(list = ls())

## ---------------------------------------
## Load Packages 
## ---------------------------------------
require('AER')
require('ivpack')
require('data.table')

## ---------------------------------------
## Load Data and Functions
## ---------------------------------------
load("../data0812.RData")

directory <- "../functions/"
functions <- list.files(directory)  
loadfunctions <- sapply(functions, FUN = function(x)source(paste0(directory, x)))

## ---------------------------------------
# Construct Instrument
## ---------------------------------------
data0812 <- constructIV(data0812)

## ---------------------------------------
#  Select Last Case Before Election
## ---------------------------------------
data0812 <- lastCase(data0812)


## ---------------------------------------
## Row 1. Main Specifications
## ---------------------------------------

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2)"
case.controls <-   "as.factor(any_drug_2) +  as.factor(any_violent_2) + as.factor(fire_arms_2) +  as.factor(any_rob_2) + as.factor(any_dui_2) + as.factor(prior_offender_2)"
demo.controls <- "age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + regis_before"

outc.1 <- "vote2012"
endo.1 <- "pti"
inst.1 <- "judgeiv"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))


m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_d <- summary(m1a1, vcov = sandwich,diagnostic = T)$diagnostics[1, 3]
m1a2_d <- summary(m1a2, vcov = sandwich,diagnostic = T)$diagnostics[1, 3]
m1a3_d <- summary(m1a3, vcov = sandwich,diagnostic = T)$diagnostics[1, 3]

row1 <- data.table(round(rbind(c(robust.se(m1a1)[c("pti"), c(1)], robust.se(m1a2)[c("pti"), c(1)], robust.se(m1a3)[c("pti"), c(1)]),
                                 c(robust.se(m1a1)[c("pti"), c(2)], robust.se(m1a2)[c("pti"), c(2)], robust.se(m1a3)[c("pti"), c(2)]),
                                 c(m1a1_d, m1a2_d, m1a3_d),
                                 c(m1a1$n, m1a2$n, m1a3$n),
                                 rep(NA, 3)
),3))

## ---------------------------------------
## Row 2. Racially Disparate Impact 
## ---------------------------------------
endo.1 <- "pti * as.factor(race)"
inst.1 <- "judgeiv * as.factor(race)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

data0812$race <- as.factor(data0812$race)
data0812 <- within(data0812, race <- relevel(race, ref = "Black"))


m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_d1 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a2_d1 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a3_d1 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]

m1a1_d2 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race)White)",3]
m1a2_d2 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race)White)",3]
m1a3_d2 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race)White)",3]


row2 <- data.table(round(rbind(# estimates
  c(robust.se(m1a1)[c("pti"), c(1)], robust.se(m1a2)[c("pti"), c(1)],  robust.se(m1a3)[c("pti"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti"), c(2)], robust.se(m1a2)[c("pti"), c(2)], robust.se(m1a3)[c("pti"), c(2)]),
  # f statistics
  c(m1a1_d1, m1a2_d1, m1a3_d1),
  # estimates
  c(robust.se(m1a1)[c( "pti:as.factor(race)White"), c(1)], 
    robust.se(m1a2)[c("pti:as.factor(race)White"), c(1)], 
    robust.se(m1a3)[c("pti:as.factor(race)White"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti:as.factor(race)White"), c(2)], 
    robust.se(m1a2)[c("pti:as.factor(race)White"), c(2)], 
    robust.se(m1a3)[c("pti:as.factor(race)White"), c(2)]),
  # f statistics
  c(m1a1_d2, m1a2_d2, m1a3_d2),
  # obs
  c(m1a1$n, m1a2$n, m1a3$n),
  rep(NA, 3)
),3))

## ---------------------------------------
## Row 3. Racially Disparate Impact (disaggregated 'White')
## ---------------------------------------

data0812$race2 <- as.factor(data0812$race2)
data0812 <- within(data0812, race2 <- relevel(race2, ref = "Black"))

demo.controls <- "age_2012 + I(age_2012^2) + Female + as.factor(race2) + vote2008 + as.factor(noteli08) + regis_before"

endo.1 <- "pti * as.factor(race2)"
inst.1 <- "judgeiv * as.factor(race2)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_d1 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a2_d1 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a3_d1 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]

m1a1_d2 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race2)WhiteA)",3]
m1a2_d2 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race2)WhiteA)",3]
m1a3_d2 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race2)WhiteA)",3]

m1a1_d3 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race2)WhiteH)",3]
m1a2_d3 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race2)WhiteH)",3]
m1a3_d3 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(race2)WhiteH)",3]


row3 <- data.table(round(rbind(# estimates
  c(robust.se(m1a1)[c("pti"), c(1)], robust.se(m1a2)[c("pti"), c(1)],  robust.se(m1a3)[c("pti"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti"), c(2)], robust.se(m1a2)[c("pti"), c(2)], robust.se(m1a3)[c("pti"), c(2)]),
  # f statistics
  c(m1a1_d1, m1a2_d1, m1a3_d1),
  # estimates
  c(robust.se(m1a1)[c( "pti:as.factor(race2)WhiteA"), c(1)], 
    robust.se(m1a2)[c("pti:as.factor(race2)WhiteA"), c(1)], 
    robust.se(m1a3)[c("pti:as.factor(race2)WhiteA"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti:as.factor(race2)WhiteA"), c(2)], 
    robust.se(m1a2)[c("pti:as.factor(race2)WhiteA"), c(2)], 
    robust.se(m1a3)[c("pti:as.factor(race2)WhiteA"), c(2)]),
  # f statistics
  c(m1a1_d2, m1a2_d2, m1a3_d2),
  # estimates
  c(robust.se(m1a1)[c( "pti:as.factor(race2)WhiteH"), c(1)], 
    robust.se(m1a2)[c("pti:as.factor(race2)WhiteH"), c(1)], 
    robust.se(m1a3)[c("pti:as.factor(race2)WhiteH"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti:as.factor(race2)WhiteH"), c(2)], 
    robust.se(m1a2)[c("pti:as.factor(race2)WhiteH"), c(2)], 
    robust.se(m1a3)[c("pti:as.factor(race2)WhiteH"), c(2)]),
  # f statistics
  c(m1a1_d3, m1a2_d3, m1a3_d3),
  # obs
  c(m1a1$n, m1a2$n, m1a3$n),
  rep(NA, 3)
),3))


## ---------------------------------------
## Row 4. Resource Deprivation
## ---------------------------------------
demo.controls <- "age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + regis_before"

endo.1 <- "pti * as.factor(quants_inc)"
inst.1 <- "judgeiv * as.factor(quants_inc)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_d1 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a2_d1 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a3_d1 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]

m1a1_d2 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)2)",3]
m1a2_d2 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)2)",3]
m1a3_d2 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)2)",3]

m1a1_d3 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)3)",3]
m1a2_d3 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)3)",3]
m1a3_d3 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)3)",3]


row4 <- data.table(round(rbind(# estimates
  c(robust.se(m1a1)[c("pti"), c(1)], robust.se(m1a2)[c("pti"), c(1)],  robust.se(m1a3)[c("pti"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti"), c(2)], robust.se(m1a2)[c("pti"), c(2)], robust.se(m1a3)[c("pti"), c(2)]),
  # f statistics
  c(m1a1_d1, m1a2_d1, m1a3_d1),
  # estimates
  c(robust.se(m1a1)[c( "pti:as.factor(quants_inc)2"), c(1)], 
    robust.se(m1a2)[c("pti:as.factor(quants_inc)2"), c(1)], 
    robust.se(m1a3)[c("pti:as.factor(quants_inc)2"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti:as.factor(quants_inc)2"), c(2)], 
    robust.se(m1a2)[c("pti:as.factor(quants_inc)2"), c(2)], 
    robust.se(m1a3)[c("pti:as.factor(quants_inc)2"), c(2)]),
  # f statistics
  c(m1a1_d2, m1a2_d2, m1a3_d2),
  # estimates
  c(robust.se(m1a1)[c( "pti:as.factor(quants_inc)3"), c(1)], 
    robust.se(m1a2)[c("pti:as.factor(quants_inc)3"), c(1)], 
    robust.se(m1a3)[c("pti:as.factor(quants_inc)3"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti:as.factor(quants_inc)3"), c(2)], 
    robust.se(m1a2)[c("pti:as.factor(quants_inc)3"), c(2)], 
    robust.se(m1a3)[c("pti:as.factor(quants_inc)3"), c(2)]),
  # f statistics
  c(m1a1_d3, m1a2_d3, m1a3_d3),
  # obs
  c(m1a1$n, m1a2$n, m1a3$n),
  rep(NA, 3)
),3))


## ---------------------------------------
## Row 5. Resource Deprivation (different baseline)
## ---------------------------------------
## Just for table purposes, change the baseline from Bottom to Mid
data0812$quants_inc <- as.factor(data0812$quants_inc)
data0812 <- within(data0812, quants_inc <- relevel(quants_inc, ref = 2))

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_d1 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a2_d1 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a3_d1 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]

m1a1_d2 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)3)",3]
m1a2_d2 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)3)",3]
m1a3_d2 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(quants_inc)3)",3]

row5 <- data.table(round(rbind(# estimates
  c(robust.se(m1a1)[c("pti"), c(1)], robust.se(m1a2)[c("pti"), c(1)],  robust.se(m1a3)[c("pti"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti"), c(2)], robust.se(m1a2)[c("pti"), c(2)], robust.se(m1a3)[c("pti"), c(2)]),
  # f statistics
  c(m1a1_d1, m1a2_d1, m1a3_d1),
  # estimates
  c(robust.se(m1a1)[c( "pti:as.factor(quants_inc)3"), c(1)], 
    robust.se(m1a2)[c("pti:as.factor(quants_inc)3"), c(1)], 
    robust.se(m1a3)[c("pti:as.factor(quants_inc)3"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti:as.factor(quants_inc)3"), c(2)], 
    robust.se(m1a2)[c("pti:as.factor(quants_inc)3"), c(2)], 
    robust.se(m1a3)[c("pti:as.factor(quants_inc)3"), c(2)]),
  # f statistics
  c(m1a1_d2, m1a2_d2, m1a3_d2),
  # obs
  c(m1a1$n, m1a2$n, m1a3$n),
  rep(NA, 3)
),3))

## ---------------------------------------
# Row 6. Prior case status
## ---------------------------------------
endo.1 <- "pti * as.factor(prior_offender_2)"
inst.1 <- "judgeiv * as.factor(prior_offender_2)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))
 
m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1a1_d1 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a2_d1 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]
m1a3_d1 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti)",3]

m1a1_d2 <- summary(m1a1, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(prior_offender_2)1)",3]
m1a2_d2 <- summary(m1a2, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(prior_offender_2)1)",3]
m1a3_d2 <- summary(m1a3, vcov = sandwich, diagnostic = TRUE)$diagnostics["Weak instruments (pti:as.factor(prior_offender_2)1)",3]

row6 <- data.table(round(rbind(# estimates
  c(robust.se(m1a1)[c("pti"), c(1)], robust.se(m1a2)[c("pti"), c(1)],  robust.se(m1a3)[c("pti"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti"), c(2)], robust.se(m1a2)[c("pti"), c(2)], robust.se(m1a3)[c("pti"), c(2)]),
  # f statistics
  c(m1a1_d1, m1a2_d1, m1a3_d1),
  # estimates
  c(robust.se(m1a1)[c( "pti:as.factor(prior_offender_2)1"), c(1)], 
    robust.se(m1a2)[c("pti:as.factor(prior_offender_2)1"), c(1)], 
    robust.se(m1a3)[c("pti:as.factor(prior_offender_2)1"), c(1)]),
  # standard errors
  c(robust.se(m1a1)[c("pti:as.factor(prior_offender_2)1"), c(2)], 
    robust.se(m1a2)[c("pti:as.factor(prior_offender_2)1"), c(2)], 
    robust.se(m1a3)[c("pti:as.factor(prior_offender_2)1"), c(2)]),
  # f statistics
  c(m1a1_d2, m1a2_d2, m1a3_d2),
  # obs
  c(m1a1$n, m1a2$n, m1a3$n),
  rep(NA, 3)
),3))

## ---------------------------------------
# Table J1
## ---------------------------------------
cont <- rbind(row1, row2, row3, row4, row5, row6)
# add column labels
colnames(cont) <- c("(1)", "(2)", "(3)")
# add row labels
cont$names <-  c("Main result: PTI",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "PTI (Black as baseline)",
                 "Std. Error",
                 "F-test", 
                 "PTI X White (Anglo & Hispanic)",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "PTI (Black as baseline)",
                 "Std. Error",
                 "F-test", 
                 "PTI × White Anglo",
                 "Std. Error",
                 "F-test", 
                 "PTI × Hispanic",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "PTI (Bottom tercile as baseline)",
                 "Std. Error",
                 "F-test", 
                 "PTI × Middle tercile",
                 "Std. Error",
                 "F-test", 
                 "PTI × Top tercile",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "PTI (Middle tercile as baseline)",
                 "Std. Error",
                 "F-test", 
                 "PTI × Top tercile",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "PTI (No prior case as baseline)",
                 "Std. Error",
                 "F-test", 
                 "PTI × Prior case",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 ""
                 )
cont <- cont[, c("names", "(1)", "(2)", "(3)")]
cat("\nPrinting Table J1: Main Effects...\n")
print(cont)


